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Abstract. This paper describes the atmospheric correction (AC) component of the Multi-Angle 
Implementation of Atmospheric Correction algorithm (MALAC) which introduces a new way to 
compute parameters of the Ross-Thick Li-Sparse (RTLS) Bi-directional reflectance distribution 
function (BRDF), spectral surface albedo and bidirectional reflectance factors (BRF) from 
satellite measurements obtained by the Moderate Resolution Imaging Spectroradiometer 
(MODIS). MAIAC uses a time series and spatial analysis for cloud detection, aerosol retrievals 
and atmospheric correction. It implements a moving window of up to 16 days of MODIS data 
gridded to 1 km resolution in a selected projection. The RTLS parameters are computed directly 
by fitting the cloud-free MODIS top of atmosphere (TO A) reflectance data stored in the 
processing queue. The RTLS retrieval is applied when the land surface is stable or changes 
slowly. In case of rapid or large magnitude change (as for instance caused by disturbance), 
MAIAC follows the MODIS operational BRDF/albedo algorithm and uses a scaling approach 
where the BRDF shape is assumed stable but its magnitude is adjusted based on the latest single 
measurement. To assess the stability of the surface, MAIAC features a change detection 
algorithm which analyzes relative change of reflectance in the Red and NIR bands during the 
accumulation period. To adjust for the reflectance variability with the sun-observer geometry and 
allow comparison among different days (view geometries), the BRFs are normalized to the fixed 
view geometry using the RTLS model. An empirical analysis of MODIS data suggests that the 
RTLS inversion remains robust when the relative change of geometry-normalized reflectance 
stays below 15%. This first of two papers introduces the algorithm, a second, companion paper 
illustrates its potential by analyzing MODIS data over a tropical rainforest and assessing errors 
and uncertainties of MAIAC compared to conventional MODIS products. 


1. Introduction 

The measured top of atmosphere (TOA) radiance is a function of both aerosol properties and 
surface bidirectional reflectance, and separating their contributions has been one of major 
obstacles for retrieval of accurate surface reflectance. The current operational MODIS aerosol 
algorithm MOD04 uses the Dark Target method (Kaufman et al., 1997; Remer et al., 2005; Levy 
et al., 2007) which relies on an empirical spectral regression coefficient (SRC) to define the 
relationship between the visible and shortwave IR (SWIR, 2.1 pm) surface reflectance. This SRC 
is optimized for vegetated and dark land surfaces where aerosols are retrieved with high 
accuracy. The lack of retrievals over bright surfaces, however, and correlation between aerosol 
optical thickness (AOT), retrieved at high spatial resolution, and surface brightness over 
heterogeneous, especially urban, surfaces (e.g., Remer et al., 2005; Lyapustin et al., 2011b), are 
known limitations of the current MODIS approach. A version of the Dark Target method has 
been implemented in the Collection 5 MODIS Atmospheric Correction (AC) algorithm MOD09 
(Vermote and Kotchenova, 2008). This algorithm makes aerosol retrievals at MODIS pixel 
resolution, and further derives surface reflectance using Lambertian assumption. Besides 
limitations of the Dark Target method, the MOD09 surface reflectance has relatively small but 
non-negligible biases from the use of the Lambertian surface model, which grow with aerosol 
loading, atmospheric mass and inverse wavelength (Wang et al., 2010; Pinty et al., 2011). 
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This paper describes the Multi-Angle implementation of Atmospheric Correction algorithm 
(MAIAC) which retrieves parameters of the bidirectional reflectance distribution function 
(BRDF) model from the time series of MODIS measurements along with spectral albedo and 
BRF without Lambertian assumption. The MAIAC cloud mask and aerosol algorithms were 
described earlier (Lyapustin et al., 2008, 2011b). Due to retrieval of SRC, MAIAC provides 
aerosol information at high spatial resolution of 1 km over both dark and moderately bright 
surfaces in a uniform manner. The AC algorithm also uses MAIAC cloud mask and column 
water vapor (Lyapustin and Wang, 2009). 

Below, section 2 gives an overview of the implemented time series processing and structure of 
the MAIAC queue storing multi-day information for the sliding window algorithm. Section 3 
describes the atmospheric correction algorithm, discusses the problem of seasonal and rapid 
surface change, and introduces the suite of MAIAC surface reflectance products. Section 4 gives 
several examples of MAIAC processing from the small- and large-scale analysis. Lastly, section 
5 gives final considerations of merits and limitations of the developed time series approach. 

Previously, MAIAC was successfully tested in applications which demand very high accuracy of 
atmospheric correction, such as Photochemical Reflectance Index (PRI) analysis (Hilker et al., 
2009. 2010). The detailed validation of MAIAC land surface products using the AERONET- 
based Surface Reflectance Validation Network (Wang et al., 2009) will be given elsewhere. A 
rather comprehensive analysis of the MAIAC land products, focused on analysis of vegetation 
greenness, and comparison with equivalent MODIS operational products for the Amazon region 
is presented in the companion paper (Hilker et al., this issue). 


2. Implementation of the Time Series Processing 

The core of MAIAC processing is a sliding window algorithm which holds between 5 (at the 
poles) and 16 (at the equator) days of imagery, depending on latitude, with up to 80 observations. 
MAIAC processing starts by gridding the received MODIS calibrated and geolocated (LIB) data 
(Wolfe et al., 1998), splitting them into tiles, and placing the new data in the processing queue 
with the previous data. Given the selected projection, the gridding process translates sensor's 
geolocated swath observations into grid cells of fixed lat-lon coordinates. In order to limit 
variation of the footprint with changing view zenith angle (VZA), the resolution of MODIS 
500m channels B1-B7 is coarsened to 1 km. The MAIAC processing uses both individual grid 
cells, also called pixels below, and fixed-size (25x25 km 2 ) areas, or blocks, required by the cloud 
mask, surface change and aerosol algorithms. To organize such processing, we developed a 
framework of C++ classes and structures (algorithm- specific Containers). The class functions are 
designed to handle processing in the various time-space scales, for example at the pixel- vs 
block-level, and for a single (last) day of measurements vs all available days in the queue, or for 
a subset of days which satisfy certain requirements (filters). 

The structure of the MAIAC processing Queue is shown schematically in Figure 1. For every 
day of observations, MODIS measurements are stored as layers for reflective bands 1-13 for the 
AC algorithm. The queue stores the retrieval results required for atmospheric correction, such as 
water vapor and aerosol information. Besides storing gridded MODIS data (tiles), the queue has 
a dedicated memory (Q-memory) which accumulates ancillary information about every block 
and every pixel of the surface for the cloud mask algorithm. It also keeps information related to 
the history of previous retrievals, for example parameters of spectral BRDF model and albedo, 
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which are required for both aerosol and AC algorithms. Given the daily rate of MODIS 
observations, the land surface represents a relatively stable background over short time intervals. 
Therefore, knowledge of the previous surface state significantly enhances the accuracy of aerosol 
retrievals and cloud detection, as well as the quality of atmospheric correction, for example, by 
imposing a requirement of consistency of the time series of BRF and albedo. 

Because MAIAC memory requirements are high, the size of the tile can be scaled to fit an 
operational memory of a particular workstation (e.g., 300-1000 km) with the constraint that it 
should be a multiple of the block size (25 km). 

3. Atmospheric Correction 

The MAIAC cloud mask algorithm has been described in detail in Lyapustin et al. (2008). 
Briefly, the algorithm composes a dynamically updated reference clear-sky image of the surface 
from spatial and time series analyses. The knowledge of reference clear-sky reflectance in 
addition to spectral and thermal reflectance tests (Ackerman et al., 1998) has been shown to 
improve cloud detection considerably (Lyapustin et al., 2008). Once the cloud mask is created 
and aerosol retrievals performed, the time series of MODIS measurements is filtered for every 
pixel and the remaining clear-sky data are placed in a “container”. The filter excludes pixels with 
clouds and cloud shadows, the snow-covered pixels as detected by the CM algorithm during 
land-water-snow classification (see Lyapustin et al., 2008), and pixels with high AOT (>1) 
where sensitivity of measurements to surface reflectance is low. The container stores 
measurements along with the corresponding radiative transfer functions from the look-up table 
(LUT). If the number of available measurements exceeds 3 for a given pixel with sufficient 
angular coverage, then the coefficients of BRDF model are computed. 

3.1 Inversion for RTLS Coefficients 

In the operational MODIS land processing, the BRDF is determined in two steps: first, the 
atmospheric correction algorithm derives surface reflectance for a given observation geometry in 
Lambertian approximation (Vermote and Kotchenova, 2008), and next, parameters of the Ross- 
Thick Li-Sparse (RTLS) model (Lucht et al., 2000) are retrieved from the time series of surface 
reflectance accumulated for a 16-day period (Schaaf et al., 2002). The retrieval is currently 
repeated once every 8 days. 

The RTLS model is also used in MAIAC. This is a linear kernel model depending on three h - 
coefficients used to describe isotropic (or Lambertian), volumetric (Roujean et al., 1992) and 
geometric-optical (Li and Strahler, 1992) terms: 

p{P 0 ,P,<P) = k L +k v f y (ii 0 ,ii,(p) + k a f a (n 0 ,n,<p). (1) 

The MAIAC AC algorithm derives RTLS coefficients directly by fitting the measured TOA 
reflectance accumulated for a period of 4-16 days. The inversion is based on equation (25) 
derived in (Lyapustin et al., 2011a): 

R( f i D ,u,<p) = R D (/i 0 ,H,(p)+k L F L (ti 0 ,!t)+k v F r (Ho,P,<P)+k a F c (n 0 ,ii,(p)+R' , ‘(n„,n). (2) 

Here, R D is atmospheric (path) reflectance and F-functions are integrals of the atmospheric path 
radiance incident on surface and atmospheric Green's function (Lyapustin and Knyazikhin, 2001) 
with respective kernels of the RTLS model. R nl is a weakly non-linear function of the surface 
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reflectance, describing multiple light scattering between the surface and atmosphere. The F- 
functions and R nl are computed analytically using eight primary functions which are stored in the 
look-up table (LUT). The fast LUT radiance-restore algorithm computes the TOA reflectance as 
a function of view geometry, wavelength, surface pressure, column water vapor and aerosol 
parameters (optical thickness and aerosol model) (Lyapustin et al., 2011a). 

Equation (2) provides an explicit parameterization of the TOA reflectance in terms of the RTLS 
BRF model parameters K={k^, k G , k y } 7 . The quasi-Iinear form of equation (2) leads to a very 
efficient iterative minimization algorithm: 

rmse = ^ (rj K) - Fj L k L(n) - Fjk V{n) - Ffk G{n) f = mjn , r (n) = R - R° - R nl{n ~ l) , (3) 

where R is measured reflectance, index j refers to different days, and n is the iteration number. 
Equation (3) provides an explicit least-squares solution for the kernel weights. In matrix form, 
the solution is written as: 

K {n) =A-'b {n \ (4) 

where 
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In the first iteration, the small non-linear term is set to zero, = 0 . and the multiple 

reflection factor a (see Lyapustin et al., 2011a) is set to one, a (0) = 1 . These parameters are 
updated once, after the RTLS coefficients are computed in the first iteration. Except for snow- 
covered surfaces, the convergence is achieved with high accuracy in two iterations in all 
conditions because the non-linear terms are small. 

Prior to inversion, the algorithm checks if the dataset has a sufficient angular sampling. This 
issue was investigated in several studies (e.g., Barnsley et al., 1997; Lucht and Lewis, 2000). The 
MODIS operational BRDF/albedo algorithm (Schaaf et al., 2002) makes an inversion if at least 7 
cloud-free observations are available during the 16-day period. Except for the dry arid regions of 
the world, this condition cannot easily be met. To increase frequency of retrievals, we studied the 
inversion problem experimentally with MODIS data varying the minimal number of 
measurements n (from 3 to 10) and testing different metrics of angular sampling. One of the 
metrics tested was the determinant of the inverse matrix A which indicates whether the sampling 
angles are sufficiently different for model inversion. Although such analysis is, perhaps, most 
straightforward theoretically, we found it often too restrictive. In the end, a simple criterion was 
chosen based on the range of cosine of the view zenith angle (AW-A^mm ^ 0.2) and a requirement 
that the measurements in both forward and backscattering range of angles are represented, which 
is usually sufficient for a robust and consistent retrieval at ri> 4. 
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3.2 Solution Selection and Update 

Although the RTLS model allows efficient inversion due to its linearity, there are several caveats 
associated with this model. First, the RTLS kernels are not orthogonal, are not positive-only 
functions, and are normalized in a somewhat arbitrary fashion that is not linked to the radiative 
transfer. These factors reduce the stability and uniqueness of the solutions, such that small 
perturbations in measurements may lead to significantly different solutions in terms of k- 
coefficients. Second, the high goodness-of-fit at the measurement angles alone does not 
guarantee the correct shape of the retrieved BRDF, and may result in negative or very large 
positive BRF values at other angles. The albedo, being an integral function, is especially 
sensitive to the BRDF shape. For these reasons, we use several tests to remove the unrealistic 
solutions. 

The initial validation of the solution checks that the maximal difference over all days between 
the measured and computed TOA reflectance does not exceed a threshold of 0.05, R Meas - R rtls 

<0.05. If it does, the day (measurement) with the highest deviation is excluded from the queue 
and the inversion is repeated. 

In the next test, the algorithm verifies that values of the direct-beam albedo ( q ) at £24=15°, 45°, 
65° are positive and that the new solution is consistent with the previous solution: 

1^(45° ) — <7 Pr "(45°)| <A(X), where A is the band-dependent threshold currently equal to 0.04 

(blue), 0.05 (green and red), 0.07 (for spectral region of 0.8 - 1 .6 fim) and 0.05 for the shortwave 
infrared band (2.1 fim). Consistency of the solution is characterized by a status index. Initially, 
the confidence in the solution is low (status^ 0). Each time the new retrieval agrees with the 
previous retrieval, status increases by 1. When status^ 3, the RTLS solution is considered 
reliable. 

The thresholds (0.05 and A(X)) in the RTLS inversion routine are selected to be strict enough to 
reject most of undetected clouds which remain the dominant source of errors, but sufficiently 
loose for the solution to adapt to potential surface change. The most pervasive type of change is 
seasonal variations, related to the spring green-up and fall senescence at northern latitudes, or 
greenness variations caused by wet and dry seasons in tropics. The total seasonal variation of 
reflectance over vegetated surface is about several absolute percent in the visible bands (-0.03- 
0.1), and is significantly higher in the near infrared (-0.1 -0.4). The threshold A(K) for the daily 
variation was selected accordingly, and our analysis of a large volume of processed MODIS data 
confirms that the MAIAC algorithm does not reject measurements when surface is changing, 
even in the agricultural regions characterized by the rapid reflectance change during green-up or 
harvesting. 

When the new solution is validated, the RTLS coefficients and direct-beam albedo #(45°), stored 
in the Q-memory, are updated. The update is done with relaxation, designed to mitigate the 
random noise of retrievals: 

K k = wK^ + (\- w )K r r . (5) 

The weight w depends on our confidence in the previous solution, which increases with its 
status. The weight w= 1 for the first retrieval (status=Q) t w=0.8 for the second retrieval, w=0.6 for 
the third consecutive retrieval, and w=0.5 thereafter. This method of update reduces random 
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noise when the surface is relatively stable, but it delays response when the surface is changing. 
To avoid the solution response delay, the status of the pixel is dropped to unity when the 
seasonal (low amplitude) change is detected, or to zero in case of detected rapid (high amplitude) 
change (see sec. 3.3). 

When the solution cannot be produced because of cloudiness or high atmospheric turbidity, we 
assume that the surface remains stable and the gaps can be filled with the previous RTLS 
solution for up to a 32-day period. During the short time intervals (1-4 weeks) when the surface 
is not expected to change, this is the most natural and accurate way of gap-filling with specific 
RTLS/albedo spectral solution for a given pixel, which does not involve additional assumptions. 
On the other hand, the gap-filling methods based on the land cover classification and other 
ancillary information (e.g.. Moody et al., 2008), may provide an alternative extrapolation for the 
longer time periods. For the user’s decision support, the gap-filled pixels are marked as 
“Extended” in the quality assurance (QA) value, with parameter QA .nDelay giving the number 
of days since the last update. 


3.3 Accommodation of Surface Change 

The AC algorithm described above works well when the rate of surface change is not high, or in 
other words, when the change in surface reflectance during the 1 6-day period does not exceed a 
certain threshold. In this case, the AC algorithm adjusts the RTLS model parameters in a timely 
fashion, as for example over the East Coast USA. On the other hand, in regions with strong 
vegetation dynamics and bright soils, e.g. USA mid- West agricultural regions, the surface albedo 
change in 1 6-day intervals may be very high. In general, the RTLS inversion during the period of 
surface change is not a valid procedure, as changes in reflectance over time are attributed to the 
particular observation angles and BRDF shape. In general, the question of how much the surface 
can change with retrievals still providing a reasonable BRDF shape, remains unanswered. For 
example, the MODIS operational BRDF/albedo algorithm adopted an empirical approach where 
the RTLS model inversion is attempted each time, and when it fails, the scaling algorithm is 
applied. In the last case, the magnitude of reflectance is adjusted (scaled) with the single last 
measurement assuming that the BRDF shape does not change significantly with the surface 
change. 

In the current MAIAC version, the change detection algorithm is triggered when the atmosphere 
is not hazy (AOTo.47<0.4 globally, and <0.6 for regions with perpetual haze, e.g. Indo-Gangetic 
plane). The change detection is based on the Red and NIR MODIS bands (1-2), and requires 
synchronous change in both bands in the opposite directions, for example decrease in Red and 
increase in NIR for green-up, and opposite change for senescence. This approach is less prone to 
errors from undetected clouds as opposed to use of the normalized difference vegetation index 
(NDVI). To remove geometry variations and enable reflectance comparison between different 
days (obtained in different view geometries), we are using the geometry-normalized BRF n : 


Pn 


= p(d 0 ,e,(p) 


R7XS(45 q ,0°) 
RTLS(d 01 e,(p) ■ 


(6) 


An example comparing BRF (see sec. 3.4) and BRF n is given in Figure 2 for a single 1 km pixel, 
covered by a forest, near Greenbelt, Maryland, USA. One can see that this approach reduces 
geometric variability by a factor of 3-6, and parameter BRF n shows strong seasonal cycles, 
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which closely resemble vegetation phenological cycles (examples for rapid surface change are 
given in the next section). 

Using the BRF n metric, we were able to answer the question of RTLS inversion during the 
period of change. Based on empirical analysis of MODIS data over different geographical 
locations, we found that the RTLS retrieval remains reasonable (does not produce BRDF shape 
artifacts, and ^-coefficients change consistently during the period of change) if relative BRF n 
change during the accumulation period does not exceed -15%, 

Ap„/p„<0.15 (7) 

According to this criterion, MAIAC has two values of change mask, ^CHANGE for gradual 
changes of 0.05<Ap w /p„ < 0.15 , and _CHANGE_BIG for abrupt, larger changes of 

A p n l p n > 0.15 . In the latter case, MAIAC uses only a scaling approach where the parameters of 

RTLS model are scaled so that the model would match the latest measurement. Such an approach 
has advantages of immediate adjustment of spectral BRDF to the surface change and faster 
computations, assuming the BRDF shape has only a secondary role. 

The detected change for individual pixels is accepted only in case of green-up which produces a 
spectrally unique signal. On the contrary, because the senescence signal can be confused with 
effect of unresolved clouds, it is only accepted if it is a large scale event covering at least a 
quarter of pixels in the 25 x25km block. 


3.4 MAIAC Surface Reflectance Products 

MAIAC computes three main land products at 1 km resolution - the set of RTLS coefficients, 
albedo and BRF in MODIS bands 1-8 and 12-13. 

The albedo is defined in the MAIAC algorithm as ratio of surface-reflected to incident radiative 
fluxes (Lyapustin et al., 201 la). Thus, albedo is derived for a given solar zenith angle in ambient 
atmospheric conditions, and can be directly compared to ground-based measurements. Given the 
^-values of the RTLS model, the so-called black-sky and white-sky albedos (Schaaf et al., 
2002), representing the hemispheric integrals from RTLS over the downward and downward and 
upward directions, respectively, under the assumption of uniform diffuse irradiance, can be 
easily calculated. MAIAC also internally computes a Normalized BRF for a fixed geometry of 
nadir view and SZA=45° (NBRF): 

NBRF>=k L k -1. 10003^ -0.04578£[. (8) 

NBRF is analogous to the MODIS NBAR (nadir BRF-adjusted reflectance) product (part of the 
MOD43 standard product suite) with an additional solar angle normalization. The spectral NBRF 
is essentially a background image of the surface useful for the product quality analysis and 
imagery applications. With the geometry variations removed, the NBRF is a useful metric for 
different applications, and users can easily compute it from Eq. (6) using provided RTLS kernel 
weights. 

The BRF by definition (e.g., see Schaepman-Strub et al., 2006) represents a reflectance which 
would be measured if the atmosphere were absent. It is produced by the MODIS operational 
atmospheric correction algorithm (MOD09) using Lambertian surface model and is commonly 
called surface reflectance. MAIAC computes BRF from the latest MODIS measurement using 


7 



the known BRDF shape. To illustrate computation of the BRF, we re-write equation for the 
measured TO A reflectance as follows: 

R(Po>P»V) m R D (9) 

where R Surf is a surface-reflected terms computed using the current RTLS parameters and 
retrieved aerosol data, and c is spectrally-dependent scaling factor. Then, 

BRF^.W) « c A R7X5 A (^ 0 ,//,<p), (10) 

where RTLS, for the direct-beam geometry is computed from the model. Because R Surf is a non- 
linear function of the surface reflectance, equation 10 is not rigorous. Yet, because the non- 
linearity is weak and scaling factor c is close to 1, the accuracy of (9-10) is high, within a few 
percent of the modelled values under most circumstances. The described scaling algorithm was 
selected in favor of an accurate computation of BRF from the direct surface-reflected term in 
order to account for variations of the MODIS footprint with the scan angle. Over heterogeneous 
surfaces, such variations can modulate the surface brightness regardless of its BRDF, and scaling 
algorithm (9-10) accounts for these variations symmetrically in the direct and diffuse surface- 
reflected terms providing overall more stable result. The BRF is computed for the land pixel if its 
RTLS parameters are known, the pixel is cloud-free and the atmosphere is not very hazy 
(A0T<1). 

A Lambertian assumption is used during the algorithm initialization period which may last from 
just 4 days in cloud-free low AOT conditions to over a month depending on cloudiness. During 
this period of time, the total algorithm performance is sub-optimal with higher rate of undetected 
clouds and reflectance biases from Lambertian assumption. In the near future, this limitation will 
be removed with automated re-run of the initialization period. 

In addition to 1 km scale, the spectral BRF in MODIS bands 1-7 is also computed at 500m 
resolution from data gridded to 500m mesh nested in 1km grid. The higher resolution BRF 
product is required in different applications such as land disturbance/recovery or change 
detection. The computation uses the same scaling approach given by equations (9-10), and 
assumes that the BRDF shape, known at 1 km scale, is also valid at 500m resolution. Obviously, 
over heterogeneous surfaces this product is much more affected by the surface spatial variability 
than the 1km product from the pure gridding prospective. This disadvantage, however, is 
outweighted by the benefits of higher resolution for variety of applications. 


4. Examples 

In this paper, only a few examples for atmospherically corrected MAI AC data are shown as the 
companion paper (Hilker et al., this issue) contains an extensive comparison to conventionally 
atmospherically corrected MODIS data over 2.88 mill, km 2 of tropical rainforest. Figure 2 
shows a time series of surface reflectance and NDVI from MODIS Aqua data for a single 1km 
pixel representing mixed forest near Greenbelt, Maryland, USA. The BRF (top) and BRF n 
(middle) data are shown for 5 spectral bands: the red (Bl), blue (B3) and green (B4) channels are 
indicated by their color, while NIR (B2) and SWIR (B7) are displayed in brown and back, 
respectively. To illustrate variability of BRF with changing view geometry, the data points for 
the NIR channel are connected by a line which shows the range of variation of 0.1-0.15 at the 
mean value of 0.35 during the maximum greenness. The BRDF normalization to the standard 
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view geometry (middle) shows a significant reduction of geometric variability, by a factor of 3-6 
for different channels. After normalization, a consistent annual pattern of reflectance becomes 
apparent even in the dark visible channels, for example a much higher seasonal variation of 
chlorophyll-related blue and red band reflectances as compared to the green reflectance. This is 
clearly visible in the BRF n time series plot (middle) while largely masked by geometric 
variations in the BRF (top) plot. The bottom plot shows that geometric normalization, at least 
during the summer season, has little effect on NDVI, where NDV1 is shown in red and NDVI n is 
shown in black. In winter, the normalization reduces the NDVI variability which is mostly 
caused by a residual snow. Several low NDVI outliers during the maximum greenness period 
shows the effect of undetected clouds, although their occurrence is generally low. 

Figure 3 gives an illustration of MAI AC performance using 50 km MODIS Aqua subsets for the 
GSFC site (Greenbelt, Maryland). This figure is a typical example of our data visualization in a 
form of time series, similarly to the data arrangement in the processing queue. It helps us analyze 
both the overall quality of solution and complex relationships between different retrieved 
parameters. Each image in Figure 3 shows 15 days of observations in 6 different columns: 1 - 
gridded MODIS Aqua TO A RGB images; 2 - MAI AC cloud mask; 3 - RGB NBRF; 4 - RGB 
BRF; 5 - MAI AC AOT at 0.466 pm with values in scale of 0-1; 6 - detected change. The blue 
and green colors in the change column indicate _CHANGE and _CHANGE_BIG greening, while the 
light blue and yellow colors stand for _CHANGE and _CHANGE_BIG browning. 

Normally, spectral BRF captures the surface change immediately while the RTLS model 
response through a regular inversion process may be delayed by several days. However, with 
detection of a big change, the implemented scaling algorithm adjusts the BRDF model 
immediately based on a single day of observations. This is illustrated in the left-hand image 
which shows results during spring green-up (DOY 92-105, 2010). While the greening change, 
indicated by a blue color, is detected from the first day, the color adjustment of the NBRF which 
represents the BRDF model, takes place on days 10-11 with detection of _CHANGE_BIG 
greening. A similar pattern can be seen in the right image showing days DOY 285-311, 2003 
during the fall senescence. Of 30 in total, the middle 15 observations are excluded as mostly 
covered by clouds. The _CHANGE_BIG browning signal, indicated by yellow, is detected in the 
last half of the shown period. 

Figure 4 shows the difference in response of the BRDF model (via NBRF) and BRF to a rapid 
surface change caused by an agricultural fire near Skukuza site, Africa. The burning takes place 
in the upper-left comer of the images with affected area marked by ovals. The burnt area rapidly 
expands during the first three days which is clearly visible in the BRF image. The change 
detection algorithm detects change, curiously, as _CHANGE_BIG greening event on day 4, with 
scaling adjustment of the BRDF model clearly visible on day 5. This example is a demonstration 
of potential of MAI AC surface reflectance products complemented with change detection mask 
for analysis of rapid surface change and disturbance (e.g., Roy et al., 2002). 

Figure 5 gives example of the large-scale MAIAC processing, showing central-eastern USA for 
September 19, 2007 at 10% of the original 1km resolution. Shown from top are the TOA RGB 
MODIS Aqua image, MAIAC cloud mask, AOTo. 466 , and finally, RGB BRF and NBRF. The 
images stitches three MODIS granules. The change of view geometry across the border of 
granules is clearly visible in both TOA RGB and BRF images. The NBRF image is computed 
from the known RTLS BRDF model for the fixed view geometry, and shows the entire gap-free 
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land surface area. The above examples show the high overall quality of MAIAC processing 
including cloud detection and correction of atmospheric effects. 


5. Concluding Remarks 

This paper is third in a series of recent MAIAC publications, describing the atmospheric 
correction algorithm, with the radiative transfer basis and aerosol algorithm reported recently 
(Lyapustin et al., 201 1 a, b). MAIAC is an attempt to build the algorithm with 
physical/mathematical model of the atmosphere-surface radiative transfer as complete as 
possible and to provide a radiatively consistent suite of atmosphere-surface properties. On one 
hand, this allows us to avoid using common assumptions, for example, about the spectral 
regression coefficient required for aerosol retrievals, and derive it from the time series and 
spatial analysis of MODIS data. The synergistic nature of MAIAC and information sharing 
among CM, water vapor, aerosol and surface reflectance algorithms benefits the overall quality 
of the product suite (see a companion paper). For example, aerosol algorithm requires spectral 
surface BRDF information, while also providing enhanced cloud detection (Lyapustin et al., 
2012). On the other hand, explicit use of the time series and a dynamic nature of the Earth 
reflectance increase complexity of the algorithm as it requires tracking surface state and 
detecting seasonal and rapid surface change for every grid cell, which we described in this paper. 

Contrary to the current MODIS operational approach, where processing is split into the AC 
algorithm (MOD09) providing Lambertian surface reflectance and BRDF/albedo algorithm 
(MOD43) retrieving coefficients of RTLS model from 16 days of Lambertian reflectances, 
MAIAC derives BRDF coefficients directly from the time series of MODIS top of atmosphere 
data, along with spectral surface albedo and BRF. 

The surface change detection in MAIAC is based on geometry normalized BRF n in the Red and 
NIR bands. The normalization to the standard view geometry of VZA=0° and VZA=45° in BRF n 
strongly suppresses geometric variations of BRF by a factor of 3-6, allowing reliable change 
detection. One of findings of this work establishes the empirical bounds for the spectral 
reflectance change during the accumulation period (ABRF n /BRF n <15%) when the RTLS 
inversion remains a generally robust procedure. When the relative reflectance change is larger, 
the AC algorithm utilizes the MOD43 scaling approach by adjusting the magnitude of BRDF 
based on a single latest measurement. 

MAIAC provides a gapless BRDF/albedo product: under cloudy conditions, it fills any given 
grid cell with the previous retrieval value for up to 32-day period assuming that the surface did 
not change, which is the most natural way of gap-filling with the pixel-specific values without 
additional assumptions. 

Because MAIAC uses a block-level processing for the SRC and aerosol retrievals (Lyapustin et 
al., 2011b), those blocks may show in the atmospherically corrected imagery at 25 km scale. 
Generally, this is a minor phenomenon with magnitude within 0.001-0.002 in reflectance units in 
the Blue-Green bands, which rapidly decreases with wavelength. However, sometimes it 
becomes apparent because our eye is trained to catch geometrically correct and reproducible 
patterns. 

Currently, MAIAC undergoes operational code conversion and testing phase. According to the 
current schedule, it will be used to process the full MODIS Terra and Aqua dataset using 
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MODIS Adaptive Processing System (MODAPS) after completion of the MODIS Collection 6 
land discipline re-processing (late fall 2012 - winter 2013). 
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Figure 1. Structure of MAIAC Queue. 

The Queue, designed for the sliding window 
algorithm, stores up to 16 days of gridded MODIS 
observations at 1 km resolution. The data are 
stored as Layers (double-indexed arrays) shown in 
the upper-left corner. A dedicated Q-memory is 
allocated to store the ancillary information for CM 
algorithm, such as a reference clear-skies image 
( refcm ), block-level statistical parameters {rl wa *; 
cri; ABT}, and results of dynamic Land- Water- 
Snow classification (mask_LWS). The Q-memory 
also stores results of previous reliable BRDF 
retrievals for MODIS bands 1-13. 
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Figure 2. A time series of surface BRF (top), BRF n (middle) and NDVI/NDVI n (bottom) from 
MODIS Aqua data for a single 1 km pixel covered by a forest near Greenbelt, MD, USA. The 
BRF/BRF n data are shown for 5 spectral bands: the red (Bl), blue (B3) and green (B4) are 
indicated by their color, while NIR (B2) and SWIR (B7) are displayed in brown and back, 
respectively. At the bottom plot, the NDVI is shown in red, and NDVI n is shown in black. 
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Figure 3. Illustration of MAIAC performance during spring green-up (left, DOY 92-105, 2010) 
and fall senescence (right, DOY 285-311, 2003) for 50 km MODIS Aqua subsets centered at 
NASA Goddard Space Flight Center, Maryland, USA. In each image, the vertical columns show: 
1 - MODIS Aqua TOA RGB images; 2 - MAIAC CM (legend: red/yellow - cloud, blue/light 
blue - clear Iand/water); 3 - RGB NBRF; 4 - RGB BRF; 5 - MAIAC AOT at 0.466 fxm with 
values in scale of 0-1; 6 - detected change (legend: pink - fill value; blue/green - 
_CHANGE/_CHANGE_BIG greening, cyan/yellow - _CHANGE/_CHANGE_BIG browning, red - 
discrepancy between the RTLS model and measured BRF in band B7 exceeds 3a). 



Figure 4. Rapid surface change caused by a fire near Skukuza, Africa. Columns and scales are 
the same as in Figure 3. The actively burning area (marked) becomes apparent immediately in 
the BRF images while the NBRF (BRDF model) response is delayed by 2 days. 
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